Leray-alpha simulations of wall-bounded turbulent flows 
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The Leray-a model reduces the range of active scales of the Navier-Stokes equations by smooth- 
ing the advective transport. Here we assess the potential of the Leray-a model in its standard 
formulation to simulate wall-bounded flows. Three flow cases are considered: plane channel flow at 
Re^ = 590, Rayleigh-Benard convection at Ra = 10^ and Pr = 1, and a side-heated vertical channel 
at Ra = 5 X 10^ and Pr = 0.7. The simulations are compared to results from a well resolved and 
coarse DNS. It is found that for all three flow cases, the variance in the velocity field increases as the 
filter width parameter a is increased, where a is connected to the filter width as = aAxi, with 
Axi the local grid size. Furthermore, the viscous and diffusive wall regions tend to thicken relative 
to the coarse DNS results as a function of a. In the cases where coarse DNS overpredicts wall 
gradients (as for Rayleigh-Benard convection and the side- heated vertical channel), the thickening 
is beneficial. However, for the plane channel fiow, coarse DNS underpredicts the wall-shear velocity, 
and increasing a only degrades the results. It is shown that buoyancy effects need to be included 
with care, because of the close relation of turbulent heat flux and the production of turbulent kinetic 
energy. 
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I. INTRODUCTION 

The usual way to surmount the computational com- 
plexity of high Reynolds number turbulent flows and 
the exorbitant grid resolution demands of high Reynolds 
number turbulent flows is to average or filter the Navier- 
Stokes equations. With such pre-processing of the equa- 
tions, one ends up with unclosed stresses which need to be 
modeled. In large-eddy simulations, the unclosed stresses 
are usually modeled with eddy-diffusivity type models. 
An alternative way of modeling is to smooth the dynam- 
ics of the Navier-Stokes equations by a direct modifica- 
tion of the advective terms [1, [H, [H, [11] . This allows 
one to derive the corresponding subscale model without 
any further assumptions, provided the employed filter C 
is invertible. The modification can be chosen such that 
important physical properties of the Navier-Stokes equa- 
tions are preserved, such as energy, helicity, the Kelvin 
circulation theorem, or symmetries. 

The Leray-a model [1, [l^, [sl] is the simplest regular- 
ization model. For this model, the advective operator 
djUjUi is replaced by djUjUi, where the filtered velocity 
Ui — C(ui) and £ is a filter operation. Historically, this 
model was used to prove existence and regularity of Ui, 
as well as convergence to the Navier-Stokes solution as 
the filter width tends to zero [1^. Analytical estimates 
invoking Prandtl-Blasius arguments Q predict an excel- 
lent match for the laminar, transitional and turbulent 



regime, although the physical relevance of the solutions 
in the turbulent regime with respect to the boundary 
conditions and turbulence intensities may be questioned 
jsot . The Leray-a model has been used for the simulation 
of the turbulent mixing layer, and showed superior accu- 
racy in comparison with dynamic eddy- viscosity models, 
in particular with respect to the quality of the small re- 
solved scales [l^. The explicit filtering concept in LES 
can be traced back to , although that model was orig- 
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inally proposed by (22|. 

The Lagrangian averaged Navier-Stokes-a (LANS-a) 
model 0, II, IS, [11,1111 in addition, poscsses 

a filtered Kelvin circulation theorem. The LANS-a 
model has been successful in simulations of homogeneous 
isotropic turbulence 0, HI] and for the temporal mixing 
layer [l3| . A regularization approach which preserves the 
skew-symmetry of the advective operator has been re- 
cently proposed by [111. Here, the energy, enstrophy (in 
2D) and helicity are conserved, and accurate results are 
obtained with coarse simulations of plane channel flow. 

Given the excellent results obtained with the Leray- 
alpha model the turbulent mixing layer [ll], it is a nat- 
ural step to assess whether the Leray-alpha model per- 
forms as well for wall-bounded flows. In this paper, we 
study three canonical flow cases which represent a wide 
range of problems occuring in practical situations. The 
first is a plane channel fiow, which is one of the standard 
test cases for turbulence models. The main production 
mechanism of turbulent kinetic energy for this flow is 
by shear production. The second test case is Rayleigh- 
Benard convection, for which the turbulence is created 
purely by buoyancy production. This is a challenging 
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test case because of the sensitive dependence of the heat- 
transfer on the turbulent dynamics. The third test case 
is a side heated vertical channel, for which turbulent ki- 
netic energy is produced both by shear and buoyancy. 
This test-case is challenging because it exhibits counter- 
gradient turbulence heat and momentum fluxes [l7l |. 



In this study, we perform a coarse DNS simulation for 
each flow case. The resolution for this simulation is cho- 
sen such that the turbulence statistics have degraded rel- 
ative to the resolved DNS solution, but the numerical 
solution is not dominated by wiggles. This simulation is 
then taken as the reference. The expectation is that the 
Leray-alpha model will improve on these results by mod- 
ifying the advective transport. The filter used to smooth 
the transport velocity is associated with the Helmholtz 
operator, as proposed in the original formulation of the 
Leray-a model [8[ . In order to avoid excessive smoothing 
in the viscous and diffusive wall regions, a non-uniform 
filter width is used in the wall-normal direction, which is 
related to the local gridsize. The non-uniform filtering 
gives rise to a nonsolenoidal filtered velocity field. As 
shown in [s^ , a non-solenoidal velocity field gives rise to 
spurious production of TKE very close to the wall, which 
degrades the performance of the Leray-a model. In order 
to circumvent this issue, a projection method is used to 
ensure that the filtered velocity be solenoidal. 



The details of the Leray-a model are outlined in section 
im A preliminary study [1^ showed that for Rayleigh- 
Benard convection, filtering with a constant filter size 
leads to an artificial thickening of the hydrodynamic 
and thermal boundary layer and significantly increased 
variances. As properly resolving viscous and diffusive 
wall region is of primary importance for wall-bounded 
flows, a space-varying filter size is applied in the wall- 
normal direction which scales with the grid resolution. 
This has implications for the velocity field, and the de- 
tails about enforcing a divergence-free smoothed veloc- 
ity field Ui when using a space-varying filter [s^l are 
discussed in section IIIII The simulations for the chan- 
nel flow, Rayleigh-Benard convection and the side-heated 
convection are discussed in section HVl The simulation of 
buoyancy-driven flows requires an extension of the mo- 
mentum equation with the buoyancy force and an addi- 
tional temperature equation. It turns out for the adopted 
formulation, the direct coupling between the turbulent 
heat-flux and the buoyancy production of turbulent ki- 
netic energy is lost, which is particularly problematic for 
Rayleigh-Benard convection (section |V| and appendix lX)) . 
Several suggestions are presented to correct for this unde- 
sired side effect. Concluding remarks are made in section 

lYIl 



II. LERAY-a REGULARIZATION 

The governing equations for the Leray-a model are 

dtUi + UjdjUi = vd'^Ui - dip + fi, (1) 
^^u, = 0, (2) 
(1 - dja'^dj)ui = Ui, (3) 

with Ui the velocity, p the pressure, v the kinematic vis- 
cosity and /; a body force. As can be seen in ([T]), the reg- 
ularization modeling approach results in a mixed formu- 
lation of filtered and unfiltered quantities. Even though 
Ui is unfiltered it should not be considered a velocity com- 
ing from direct numerical simulation (DNS) because of 
the modified nonlinearity in ([T]). Therefore, both u, and 
Ui are regularized variables. The variable a^ corresponds 
to the filter width which can have different values per 
direction and varies in the wall-normal direction. The 
filter C is given by the inverse of the elliptic equation 
and boundary conditions; we will denote the filtering 
operation as 

u, = C{u,) ^ {1 - djapj)-' u,. (4) 

This is the standard filter used for analytical studies of 
the Leray-a and LANS-a model [1, [l^ . Naturally, other 
filters could be used as well, as the effect of the filtering 
would (hopefully) be restricted to the small scales only. 
One advantage of using a formulation like ^ is in the 
natural treatment of boundaries, as explicit filters require 
special treatment near the boundary. However, solving 
([3]) is quite costly, and care should be taken to ensure the 
incompressibility of the filtered velocity u [1^ . 

The Leray-a model does not add additional dissipation 
to the equations. Indeed, the equation can be written in 
an LES template as 

dtUi + djUjUi + diP - lydjUi = d^ruij, (5) 

where rriij is the subfilter model, which is given by 
niij = u'jUi. Here, = Uj — Uj is the high wavenum- 
ber part of Ui. Multiplying with Ui results in 
UidjiTiij = dj{uj^uf), which shows that energy is only 
redistributed and not dissipated. 

One can also derive the implied subfilter scale model 
of the Leray-a model in terms of the filtered velocity u. 
Filtering ([!]) with £ leads to an explicitly filtered LES 
template, given by 

dfUi + djUjUi + dip — vd^Ui ~ djfhij, (6) 

but with a non-standard asymmetric subfilter term 

fhij = UjUi — UjUi ~ UjUi — UjUi + UjUi (7) 

The Leray-a model reduces the number of active scales 
by slowing down the energy cascade at small scales. This 
is shown schematically in Fig. [1] for the energy spectrum 
E of unfiltered velocity. The wavenumber associated with 
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FIG. 1; Sketch of a typical spectrum for the unfihered velocity 
Ui of the Leray-a model. 



the filter size is ka- For k < ka, the energy cascade 
is unaffected which results in the classical E oc k~^^^ 
spectrum. However, for k > ka, the filter affects the 
cascade time-scale and the spectrum is modified to E (x 
Because of the higher variances at small scales, 
the dissipation will take place on larger scales than the 
Kolmogorov scale with associated wave number fc^ . Note 
that the energy density spectrum of the filtered velocity 
E = C'^{k)E(k) falls off steeply as oc fc^"/3 beyond 
ka [see alsoQ. 



III. ENFORCING A NONDIVERGENT 5, FIELD 

Enforcing the incompressibility condition djUj ~ 
does not automatically ensure that djUj = for a wall- 
bounded flow. For uniform a, the solution /(x) = djUi of 
the homogeneous differential equation f — djOjdjf = is 
nontrivial and does not vanish upon applying the bound- 
ary conditions [s^l. The main problem with djUj ^ 
is that the purely redistributive character of advcction is 
lost, which can be demonstrated by considering the effect 
on the balance of the turbulent kinetic energy (TKE). 
The turbulent transport of TKE by fluctuations can be 
calculated by taking the advcction term djiu'^u'^) from 
the equation of velocity fluctuations , multiplying it by 
Ui and averaging: 



Kclj{uWi) 



(8) 



with e' = ^u^M^. The first term on the right-hand side is 
in divergence- form and hence purely redistributive. The 
second term is a production/destruction term, which nor- 
mally vanishes because the fluctuating field is divergence 
free. However, when djUj ^ this term can become 
nonzero, thereby allowing for spurious TKE production 
and destruction. As was shown for Raylcigh-Bcnard con- 
vection at Ra = 10^ and Pr = 1, the term e'dju'j is 
responsible for significant energy production in the hy- 
drodynamic boundar y la yer (up to 25% of the volume- 
averaged production) [33 • 

A nonuniform filter width a, creates another oppor- 
tunity for djUj ^ 0. Applying the filter C to the di- 



vergence of the unfiltered velocity gives that C{djUj) = 
djC{uj)-\-[C, dj\uj, where the brackets represent the com- 

the 



mutator defined as [A, B] = AB — BA. Therefore, 
divergence of the filtered velocity is given by 



(9) 



This shows both ways in which u can cease to be 
divergence-free. First the commutator may not vanish 
because of a nonuniform filter width a(x). Second, the 
term £(0) is nonzero when the solution /(x) = djUi of 
the homogeneous differential equation / — djOjdjf = 
does not yield / = after applying the boundary condi- 
tions, as discussed above. In |32l |. several methods were 
proposed to maintain a nondivergent Ui field: enforcing 
djUj = directly (and thus dropping djUj = 0), letting 
a vanish at the wall, using free-slip boundary conditions 
for Ui and introducing a projection method. Of these 
four options, only the projection method can be success- 
fully applied for nonuniform filters q;(x), so this is a con- 
venient choice. The projection method [e.g. [ll| adds a 
gradient dicj) to (g]) as Ui+dtcf) = C{ui) and solves a Pois- 
son equation 9^0 = diC{ui) to make Ui incompressible. 
The boundary conditions imposed on Ui and are equiv- 
alent to the conditions imposed on and p; Ui ~ and 
dn(t> — 0, respectively. 

With the projection included, the filtered velocity Ui 
is related to the unfiltered velocity Ui by a filter Q as 



g{u) = o c{uj), 



(10) 



where Pij is the projection operator. Note that Q takes 
the complete velocity vector as its argument, as opposed 
to £ which only requires one velocity component at a 
time. 

The grid resolution is normally a good indication for 
the location of the most demanding flow features. There- 
fore, we couple the flltcr size ai directly to the grid res- 
olution by a parameter a as 



(11) 



Note that the filter size is anisotropic. In the wall-normal 
direction, the grid is clustered towards to wall so that the 
filter is small near the wall and larger in the midplane. 
The smoothing parameter a is varied between < a < 1, 
for which the typical filter width spans up to four grid 
cells. This can be shown by calculating the filter width 
A from [11: 



1 

A 



Glix,x')dx', 



(12) 



where Ga is a normalized filter kernel. The free-space 
Green's function associated with the Helmholz equation 
is 



Ga{x,x ) = —exp I 

2a V a 



(13) 
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Substitution of the above expression and (fTTj) into (fT2|) 
leads to A ~ 4a Ax^. 

The projection method has been implemented in our 
code for direct simulation as a two-step process. In the 
first step, the elliptic equation ^ is solved by a direct 
method which takes advantage of the homogeneous direc- 
tions. In the second step, the projection method is used 
to make this field divergence free. The code is based on 
staggered finite differences and uses second order central 
differences in space and a second order Adams-Bashforth 
scheme in time (soj . The code is fully parallellized and 
supports grid clustering in wall-normal direction. Special 
care has been taken to preserve the purely redistributive 
character of the advection scheme on the nonuniform grid 
by using a symmetry-preserving discretization [s^ . 

IV. RESULTS 
A. Plane channel flow 

The first wall-bounded flow case under consideration 
is plane channel flow. Here we compare results for the 
Leray-a model with direct simulations [131 at Re,- = 590. 
The coordinate system is defined in Fig. [5^, with the 
X— , y— and z-coordinate aligned with the streamwise, 
wallnormal and spanwise direction, respectively. The 
medium is air with a kinematic viscosity v = 10~^m^/s 
and the domain size x Ly x Lz is 2TrS x 26 x tt6 where 
(5 = 1 m is the channel halfwidth. No-slip boundary con- 
ditions are applied to the top and bottom plates, and 
periodic boundary conditions are used on the sidewalls. 
The average velocity is fixed at U = 0.115 m/s for all 
simulations. At fuU DNS resolution (384 x 384 x 256) 
cells, this flowrate results in a shear Reynolds number 
Rct- = 610. For the Leray-a and coarse DNS simula- 
tions, the grid has been chosen such that the grid is too 
coarse for a good solution without a subgrid model, but 
fine enough so that the solution is not dominated by wig- 
gles. For the plane channel flow this resulted in a grid 
of 64 X 64 X 32 cells with strong grid stretching (up to 
16 percent near the walls). All results were averaged for 
about 60 typical turnovers (based on a typical turnover 
time t* = U/{25). 

Three different values of the filter parameter a are used 
here: a = 0.25, 0.50 and 1.00. As discussed before, 
the filter is proportional to the local grid resolution as 
Ui — aAi'i, and due to the grid clustering in the wall- 
normal direction, the filter is much smaller near the wall 
than in the midplane (Fig. [3]). In this way, the artificial 
thickening of the viscous wall region [viscous sublayer 
and buffer layer |2^. due to the filtering which hampered 
preliminary simulations [sst is minimized. 

Shown in Fig. ^Bp, are the average velocity profiles in 
plus-units, defined as u+ = u/ut and = yur/v- Here 
Ut is the friction velocity defined as = t^/ p, where 
Tw = P'dyu\^ is the shear stress at the wall. As can be 
expected, the coarse DNS velocity profile lies above the 



reference profile, which reflects an underestimated wall- 
shear stress. The results for the Leray-a model are de- 
noted by the circles, squares and crosses for a = 0.25, 
0.50 and 1.00 respectively. The simulation with a = 0.25 
has little influence on the average velocity, but as a in- 
creases the profiles deviate more and more from their de- 
sired value. The source of these deviations is a decreasing 
wall-shear stress. Shown in Fig.[2t is the shear-Reynolds 
number Rcr = u-rbjv as a function of a. The dashed 
line corresponds to the expected value (Re,- = 590). The 
coarse DNS results (a = 0) show that the insufficient res- 
olution results in a 10% underestimation of Re^. The 
subfilter scale model would ideally compensate this ef- 
fect. However, as a is increased, Re,- decreases up to an 
underestimation of 30% at a = 1.00. 

The variance profile of the streamwise velocity compo- 
nent, u'm', is shown in Fig. [^Ji, nondimensionalized by 
M^. The coarse DNS overpredicts the peak of u'u'^ by 
a factor two, and increasing a only makes the situation 
worse, with a factor five overprediction at a = 1.00. It 
should be noted that part of this increase is due to the 
normalisation: the decrease in the friction velocity causes 
u'u'^ to become larger. However, even without scaling, 
the variance increases a factor 2.5 when comparing the 
resolved DNS results to the Leray-alpha simulation at 
a = 1. 

The peak in u'u'^ can be seen to shift outwards, indi- 
cating a thickening of the viscous wall region. The proba- 
ble mechanism for a thickening of the viscous wall region 
resides in the modification of the momentum-flux v'u' in 
the buffer layer. The filtering of v will damp near-wall 
fluctuations, thereby reducing the turbulent momentum 
flux. As a result, the turbulence transport takes over 
further away from the wall than without filtering. 

To get an idea of the typical lengthscales that cause 
the overprediction in the variance, the spatial spec- 
trum of u'u' is presented. The spectrum is obtained by 
one-dimensional Fast Fourier Transforms (FFT) in the 
streamwise (spanwise) direction, and averaging over the 
other homogeneous direction. In addition, the spectrum 
is averaged over about 60 typical timescales to eliminate 
slow transients. In Fig. [^Is, the spectrum of the stream- 
wise velocity u is shown in the midplane of the chan- 
nel. The coarse DNS slightly underpredicts the variance 
on the large scales and overpredicts the variance at the 
intermediate wavenumbers. When a is increased, first 
the variance increases at intermediate wavenumbers. For 
larger a, a significant increase in the variance occurs for 
the large wavenumbers as well. It should be noted that 
in the midplane, the difference between the DNS and the 
Leray-a model is relatively small; in the viscous wall re- 
gion, the difference would be even greater. 

More insight into the enhanced variances may be ob- 
tained by studying the budget for turbulent kinetic en- 
ergy (TKE). The equation of TKE can be obtained by 
multiplying the fluctuating part of ([1]) by u[ and averag- 
ing over the homogeneous directions and over time. The 
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FIG. 2: The effect of varying the filter size a for plane channel flow, (a) definition sketch; (b) average velocity; (c) shear- 
Reynolds number Re^; (d) streamwise variance u'u'; (e) streamwise spectrum Euu in midplane; (f) production of turbulent 
kinetic energy. 



average • is denoted by an overlinc. This results in 



= —v'u'dyU — v{dju'j)[dju'^) — dy{v'e' — udyC + v'p'), 

V e T 

(14) 



where e = \u'jU[ and e' = \u[u[] V ^ e and T repre- 
sent production, dissipation and transport of turbulent 
kinetic energy respectively. The smoothed transport ve- 
locity directly modifies the shear production term V and 
the transport of the velocity fluctuations v'e' . 
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FIG. 3; The distribution of a as a function of y. 

The fact that the variance increases while the produc- 
tion of TKE remains roughly the same is one of the in- 
triguing features of the Leray-alpha model. Because of 
the attenuation of small-scale dynamics (Fig. [T|), it is 
to be expected that the total variance will increase due 
to a slowing down of the cascade at high wavenumbers. 
A recent study suggests that the variance accumula- 
tion at subfilter scales is even higher; for the LANS-alpha 
model, the small scales seem to behave as '"rigid rota- 
tors'" which are advected passively by the larger scales. 

The production of TKE is shown in Fig. [2t, nondimen- 
sionalized by / v. The thickening of the viscous wall re- 
gion is clear in the profile of , where the peak (which 
occurs where the viscous stress is equal to the Reynolds 
stress) shifts from about y+ — 15 for the coarse DNS to 
around j/+ = 30 for a = 1.00. Note also that the width 
of the production peak increases: only at y+ = 80 is P 
at the value of the reference DNS for a = 1.00. 



B. Rayleigh-Benard convection 

Rayleigh-Benard convection (RB) is generated when a 
fluid in between two flat plates is heated from the bottom 
and cooled from the top (Fig. 3^). The system can be 
characterized by the Prandtl number Pr = i/n'^^ and the 
Rayleigh number Ra = /3gAQH^{i'K)^^ . The system re- 
acts by a convective motion which is characterized by the 
Reynolds number Re = UHv~^ and by an enhanced heat 
transfer through the Nusselt number Nu = (f>H {k/S.Q)~^ . 
Here U \s& characteristic velocity and the realised heat- 
flux at the wall. Both Re and Nu are non-trivial func- 
tions of Ra and Pr and are still the subject of ongoing 
research [e.g. [l|. The coordinate system is defined with 
the z-direction pointing upwards and the gravity vector 
is in the negative z-direction. 

The medium has most material properties of water ex- 
cept for the Prandtl number, which is taken as Pr = 1 
instead of 7 to relax the computational demands [s^ : 



the viscosity v = 1.07 x 10 ^ m^/s, expansion coefficient 
/3 = 1.74 X 10-** A'-i, Ae = 2 if, with Tq = Ae/2 
and Ti = — A0/2. The domain size x Ly x is 
TH X TH X H, with = 0.15 m and the aspect ratio 
r = 4. Fixed temperature and no-slip velocity boundary 
conditions are enforced at the top and bottom plates. 
Periodic boundary conditions are used for the sidewalls. 

The simulation of a buoyancy driven flow requires 
the extension of the Leray-alpha model with a trans- 
port equation for temperature. This raises the question 
whether the advection of scalars should be with the fil- 
tered or the unfiltered velocity. Here the same modi- 
fied advection operator will be used for all transported 
quantities, i.e. advective transport with Ui. Invoking the 
Boussinesq approximation, the effect of buoyancy can be 
included by introducing a body force (3giQ in the momen- 
tum equations only. After the extension with a transport 
equation of temperature and with the additional body 
force, the Leray-a model becomes 

dtUi + UjdjUi = lyd^Ui - dip - f3giQ, (15) 
dte + UjdjO = ud^e (16) 

In the current coordinate system, the gravity is in the 
negative z-direction as gi — —gSis. Due to the homogene- 
ity and the absence of a forcing in the x and y direction, 
Ui — 0, and the system is statistically one-dimensional. 
One of the characteristic features of RB is that the total 
vertical heat-flux through the fluid is constant. Aver- 
aging over the homogeneous directions, integrating (jl6[) 
with respect to z and substituting the Dirichlet boundary 
conditions for temperature results in 

iFG'- K9^e = ^^TVw. (17) 
H 

Note that the turbulent heat-flux is in terms of the fil- 
tered velocity Ui. However, the TKE production by buoy- 
ancy is given by w'Q' (|18p and this mismatch will prove 
to be of importance, as discussed below. 

Although many results can be found in the literature 
on direct simulation of Rayleigh-Benard convection in 
wide aspect ratios [13, [2l|j no database is available as 
for plane channel flow and side-heated convection simu- 
lations. Therefore, the results will be compared to our 
own DNS results [3l|, [H at Ra = lO'' and Pr = 1. The 
grid for the DNS is 256 x 256 x 256 cells which is of suf- 
flcient resolution not to require grid clustering near the 
wall. The grid resolution for the Leray-a and coarse DNS 
simulations is 80 x 80 x 64, which is again chosen such 
that it is too coarse for DNS but is not yet dominated by 
numerical contamination. Here, grid clustering is applied 
such that of the 64 wallnorm al p oints, 8 are present in 
each thermal boundary layer [41| 

One of the most important integral flow properties in 
RB is the Nusselt number Nu. At the current Ra and 
Pr, the DNS gives that Nu = 16.1. The coarse DNS 
significantly overpredicts this value with Nu = 20 (Fig. 
|3|d), due to the insufficient resolution. The infiuence of 
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FIG. 4: The effect of varying the fiher size a for Rayleigh-Benard convection, (a) definition sketch; (b) Nusselt number Rer; 
(c) wall-normal variance w'w'; (d) wallnormal spectrum E^w in midplane; (f) production of turbulent kinetic energy. 



the Leray-a model is to decrease Nu, and at a = 0.50 
the Nusselt number Nu is approximately at its expected 
value. 

In Fig. Ilh, the wall-normal velocity variance w'w' is 
shown, normalized by f/^ where U = ^/JJgAQH is the 
free-fall velocity. As for the plane channel flow, the vari- 



ance can be seen to increase as a becomes larger. Spatial 
spectra are collected by performing a 2D FFT and in- 
tegrating over circles fc^ + ky = . Shown in Fig. [Hi 
is the spectrum of the vertical velocity in the midplane. 
The increase of variance for the coarse DNS seems to be 
mainly concentrated at the large scales (low wavenum- 
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bers). When a is increased, it can be seen that the vari- 
ance at the intermediate wavenumbers increases. At the 
low wavenumbers, the trend is in the right direction. 

As before, wc study the equation of TKE, which is for 
Raylcigh-Bcnard convection given by 

= f3gu7W - lAdJ^JidJt^) - d,{w^ - i^d^e + IJ/j/) . 

V e T 

(18) 

For this flow case, there is no production on average of 
turbulent kinetic energy by shear and the only effect of 
the filtering is in a modified transport of the velocity 
fluctuations w'e' . 

The production term V is shown in Fig.|4j;, nondimcn- 
sionalized by {PgAQHf/^/H. The coarse DNS overes- 
timates V ^ which is consistent with the overestimation of 
Nu as these are directly coupled for RB by the exact rela- 
tion [see e.g.[ll|2| {V)^ = (e)^ = ^Ra{Nu - l)Pr-^. 
The total production decreases for a = 0.25 and a — 0.50, 
but for a — 1.00, V increases again. This is not con- 
sistent with the monotonically decreasing trend for Nu 
(Fig. |3|d). Furthermore, the shape of the production 
profile changes: the constant production in the bulk is 
replaced by a production profile which peaks near the hy- 
drodynamic boundary layers. Both the trendbreak and 
the change of shape arc the results of a disparity between 
the buoyancy production term V = (3gw'Q' and the tur- 
bulent hcat-fiux w'Q' . This is a fundamental issue which 
will be discussed in detail in section fVl 



C. Side- heated vertical channel 

The side-heated vertical channel is a case where both 
buoyancy and shear production are important. The fiow 
has several unusual features, such as negative shear- 
production in the boundary layers [i^ and counter- 
gradient heat fiuxes. The side-heated vertical channel 
has been studied intensively both experimentally 0] and 
by direct simulation [3, 113, [H, |3§| . The flow geometry 
for the side-heated vertical channel is sketched in Fig. [5^. 
For this case, the x-dircction in the wall-normal direction 
and the z-direction is pointing upward. The Leray-a sim- 
ulations will be compared to the DNS database of [1^. 
The size of the domain is x Ly x = H x 6H x 12H. 
The unusually large domain size is required to prevent 
long-range correlations from influencing the statistics 
[3^ . The medium is air with a viscosity v ~ 1.0 x 
10~^ m^/s , expansion coefficient /3 = 3.3 x 10^^ and 
a Prandtl number Pr = 0.709. The distance H between 
the plates is H ^ 0.2 m and with a temperature differ- 
ence Ae = 2.7 K with Tq = Ae/2 and Ti = -Ae/2. 
The Rayleigh number Ra = (3g/^<dH^{vK)^^ for this case 
is Ra = 5 X 10^. The resolution for the coarse DNS is 
64 X 96 X 192, which as before is chosen such that the 
grid is too coarse for accurate predictions but flne enough 
to prevent that the results are dominated by numerical 
contamination. The statistics have been collected over 



25 typical turnover times. 

In accordance with [33, HI], a body force is intro- 
duced which ensures a zero mass-flux. The advantage 
is that this suppresses slow transients, thereby reducing 
the required simulation time. In addition, for high a, 
the Leray-a model without body force the mass-flux be- 
comes nonzero, which obfuscates comparison with other 
simulations. 

The average velocity proflle is shown in Fig. [5)o, nor- 
malized by the free-fall velocity ^/(3gAQH. The coarse 
DNS overestimates the average flow velocity in the chan- 
nel. When increasing a, the velocity decreases and 
the flow profile is approximately correctly predicted for 
a = 0.25. For a = 0.50, the velocity is underpredicted. 
At a = 1.00, the velocity profile remains at roughly the 
same amplitude as for a = 0.5 but here a 'kink' can be ob- 
served in the mean velocity around x/H ~ 0.2, which is 
not present in the other simulations. This kink is respon- 
sible for the shift of the shear production term u'w'dxW 
from the center to the near-wall region. 

The results for Nu = dxQ\wH / kIS.Q as a function of a 
are shown in Fig.[5j3. The increase in the average veloc- 
ity does not infiuence Nu very much for the coarse DNS. 
The trend for increasing a is that Nu decreases, similar 
to Rayleigh-Benard convection. Around a = 0.25, Nu is 
at its expected value, and above this value, Nu is under- 
predicted. The shear Reynolds number Re,- is defined as 
Rcr = UtH/v. The infiuence of a is relatively small for 
this flowcase (Fig. [S}!), with an undcrprediction of Rc,- 
of 10% at a = 1.00. 

The streamwise velocity variance proflle w'w' is shown 
in Fig. [St, nondimcnsionalizcd by f3gA&H . The coarse 
DNS overestimates the variance by 50% in the bulk. Here 
the effect of increasing a is to decrease the variance for 
a = 0.25. As a is increased further, the variance becomes 
practically constant in the bulk (a = 0.50), after which 
the maximum shifts from the center to the near-wall re- 
gion (a 1.00). 

The equation of TKE for the side-heated vertical chan- 
nel is given by 



— Pgw'Q' — u'w'dxW — v{dju[)[djU^) — dx{u'e' — vdxS. + u'p') . 

V e T 

(19) 

For this flow case, there is production of TKE both by 
shear and by buoyancy. Both components of V are shown 
as a function of z/H in Fig. [5f. Here, it can be seen that 
the lack of resolution mainly affects the shear-production 
term —u'w'dxW through the overestimation of the aver- 
age velocity (Fig.[5j3). The effect of increasing a reduces 
the shear-production and the production is estimated ap- 
propriately at a = 0.50. At a = 1.00, the buoyancy pro- 
duction term changes little with respect to a = 0.50, but 
the shear production term changes dramatically. 
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FIG. 5: The effect of varying the fiher size a for side-heated convection, (a) definition sketch; (b) average velocity; (c) Nusselt 
number Nu; (d) shear- Reynolds number Re^; (e) streamwise variance w'w'; (f) production of turbulent kinetic energy. 
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V. DISCUSSION 

The simulations show two major trends. First, gra- 
dients at the wall (as reflected in integral quantities like 
Rbt and Nu) tend to decrease. Both for Rayleigh-Benard 
convection and the side-heated vertical channel, coarse 
DNS overpredicts wall gradients, and the Leray-a model 
can improve results. In the case that coarse DNS under- 
predicts gradients at the wall, as for the plane channel 
flow, the Leray-a model does not improve the results. 

Second, the variances tend to increase as a function of 
a, specifically at the low and intermediate wavenumbers. 
In view of the model spectrum (Fig.[T]), which predicts an 
increase in the variance at wavenumbers larger than ka 
only, this may be quite surprising. However, this seems 
to be an intrinsic property of the Leray-a model. Simula- 
tions with the Leray-a model in absence of walls [l^, [l3| 
also results in enhanced variances and the low and in- 
termediate wave lengths. The presence of walls does not 
change this property. However, the presence of a wall 
seems to enhance this feature of the model, due to an 
increase in the turbulent shear production term. 

In simulations with coarse grids, a significant part of 
the fluctuations are sub-grid, and separate modeling may 
be required [e.g. Therefore, a pragmatic solution to 
remedy the increased variances may be to add some extra 
diffusion to the model. This may be done in several ways: 
1) by using a simple eddy- viscosity model; 2) by using a 
dynamic Smagorinsky procedure; and 3) by using a spec- 
tral dissipation procedure. It is clear that this procedure 
would have to be done with care, as one could easily over- 
whelm the Leray-a contribution to the subfilter stress. 
The tensor-diffusivity model (also called Clark model) 
[40| . is a good example of a model which significantly 
benefits from some additional dissipation. This filter re- 
produces approximately 90% of the subfilter stresses in 
a priori studies . However, in a posteriori studies, the 
tensor diffusivity model requires additional dissipation 
[40| . Furthermore, extra dissipation may be unavoidable 
to be able to use Leray-a with relatively coarse meshes. 
In this study, the resolution was chosen such that coarse 
DNS simulations were incapable of reproducing correct 
statistics but were not dominated by wiggles. When even 
coarser grids would have been used, the absence of addi- 
tional dissipation leads to a deterioration of the results 
of the Leray-a model. 

In this study we chose to keep the filter as close as 
possible to the original formulation The filtering re- 
quires inverting a Helmholtz equation for all three veloc- 
ity components, which is computationally expensive. In 
addition, the projection step to make divergence-free 
requires solving another Poisson equation. Even though 
the computational overhead for solving a Poisson equa- 
tion is minimized by using FFTs in the homogeneous 
directions, this procedure is more expensive than conven- 
tional SGS models. To reduce the number of operations, 
one could use explicit filters [3, l23 or apply ^ but with 
only one or two Jacobi iterations [35|. 



An important question for future research is related 
to the influence of the filter type. It is not inconceiv- 
able that an explicit filter (i.e. a filter that docs not re- 
quire solving a Poisson equation, such as a top-hat filter) 
might perform better than the Helmholtz filter. Indeed, 
the Helmholtz filter suffers from C{Q) ^ even for uni- 
form a, which violates the filtering framework (3^ . A 
projection method was required to remedy this problem. 
In addition, the Helmholtz filter does not have compact 
support, causing the filtered velocity to be influenced rel- 
atively strongly by large fluctuations far away. 

The inclusion of the buoyancy term in the Leray-a 
model has to be done with care. In the direct simula- 
tions of Rayleigh-Benard convection, the average buoy- 
ancy production V = /Sgw'Q' is constant in the bulk, 
as V depends directly on the turbulent heat flux w'Q' 
(which is constant in the bulk). This is not the case for 
the Leray-a model, where a peak in V is created near the 
hydrodynamic boundary layer (Fig. HJi). This change is 
a direct consequence of the modification of the turbulent 
heat flux. In the present formulation (|15m6p . the direct 
coupling between TKE production and heat-flux is bro- 
ken, as the equations for TKE and average heat-flux are 
given by (see also appendix TK^ : 



dz w'e' + p'w' -ud^e] ^ (igw'W - e, (20) 



w'Q' - Kd^e = ——Nu. (21) 



The buoyancy production term is given by (3gw'Q' , while 
the turbulent heat flux (which is constant in the bulk) is 
given by w'Q' (fT7|) . The difference between w'Q' and 
•w'W can be calculated by substituting w ~ w — djajdjW 
(and therefore not correcting for compressibility effects 
discussed in Sec. IIIip into w'Q' results in 

= - (ale'd^w') + a^ld~W)(dJw^. (22) 

Interestingly, the two extra terms on the right-hand side 
correspond to terms typically encountered in the trans- 
port equation of w'Q' The first is normally associ- 
ated with the molecular diffusive transport, and the sec- 
ond with the dissipative cross-correlation term, although 
the sign is opposite here. 

For buoyancy driven flows, the variation of a as a func- 
tion of the wall-normal coordinate z seems be of impor- 
tance for the variation of V over the vertical in the bulk. 
Indeed, the preliminary simulations with constant a (and 
free-slip conditions for u) did not have a variation of V 
over the height [33|. At a = 0.5, the difference of V be- 
tween the peak value at the edge of the boundary layer 
and the core is 10% (Fig. UJ;). This is quite a large dif- 
ference, and it seems worthwhile to explore whether this 
effect can be circumvented. The first and most obvious 
way is to experiment with different filters. Another op- 
tion is to modify the temperature equation. Instead of 
holding on to the same advection operator for all trans- 
ported quantities, one could apply the modified advection 
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operator Ujdj to the momentum equations only. For the 
other transported quantities, Ujdj could be used, with 
the understanding that Uj is also a regularized veloc- 
ity. This would directly restore the coupling between the 
turbulent heat-flux and the buoyancy production. How- 
ever, the danger in a formulation like this may be an 
excess of variance at small scales, which results in a high- 
wavenumber forcing contaminating the simulation. 
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VI. CONCLUDING REMARKS 

Numerical simulations of the Lcray-a model have been 
carried out for plane channel flow, Rayleigh-Benard con- 
vection and the side-heated vertical channel. The sim- 
ulations have been compared to DNS and coarse DNS. 
In general, the simulations show two trends. First, the 
viscous (and diffusive) wall region tends to thicken as a 
function of the filter width parameter a, which causes 
wall gradients such as the shear stress or the heat-flux 
to decrease. When coarse DNS overpredicts wall gradi- 
ents, such as for Rayleigh-Benard convection and (to a 
lesser extent) the side-heated vertical channel, the Leray- 
a model can improve the results. However, when the gra- 
dients are initially underestimated, such as for the plane 
channel, results do not improve. Here, additional wall- 
modeling may be needed to enhance turbulence levels in 
the near-wall region. Second, the variance at low and 
intermediate wavenumbers increases upon increasing the 
fllter size parameter a. This leads to overpredicted vari- 
ance in the velocity field which is undesired. 

An important point for buoyancy driven fiows is how 
to include a temperature forcing into the Leray-alpha 
model. Indeed, the intuitive extension (|15m6p causes the 
production of TKE by buoyancy to be no longer directly 
coupled to the turbulent heat-flux. For Rayleigh-Benard 
convection, this leads to a varying TKE production in 
the bulk (and as a function of a). 

In this paper, the performance of the Leray-a model 
was assessed for three wall-bounded flows. The Leray-a 
model was implemented in its original form, i.e. with 
the Helmholtz filter. This study indicates that, within 
this formulation, the potential of the Leray-alpha model 
is rather limited for wall-bounded flows. Indeed, the 
overpredicted variances, in particular the accumulation 
of energy on the low and intermediate wavenumbers 
pose a challenge for accurate predictions with the Leray- 
a model. A study on alternative filter types, poten- 
tially remedying the aforementioned downsides of the 
Helmholtz filter, would be a valuable next step. 
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Thc relations relating the volume-averaged dissipation 
rate of kinetic energy (e)^ and temperature variance 
(ee)^ to the Rayleigh number Ra, the Prandtl number 
Pr and the Nusselt number Nu are given by [il,!!^ 



(£), = -^Ra(Nu-l)Pr-2, 



(Al) 
(A2) 



where (e)^ and and (£e)z = (ee^av)^ + (ee^fl)^ are defined 
as 



(e)^ = z.^(aX)(a,<)^ 

(£e,av), = «;((9.e)2),, 

Here, (•)^ represents the averaging over the entire fluid 
layer, and the combination {') ^ is equal to a volume- or 
ensemble-average, when the system is ergodic. Below we 
will derive the relations for (s) ^ and (se) ^ for the Leray-a 
model. 

The equations (jAl[) and (jA2[) are obtained from the 
equations of the vertical heatflux, turbulent kinetic en- 
ergy and the two equations for temperature variance. In 
the case of the Leray-a model, these equations are given 

by 



H 

dz{w'e' + p'w' — vdzE^ = jSgw'Q' - 



(A3) 
(A4) 



dz ( w'Q' e - Kdz-Q e ) = w'WdzQ - ee,av, (A5) 



(A6) 



4< 



Averaging these expressions over the height, we obtain 



lii'Q') ^^(Nu-1) 



(A7) 
(A8) 



ft2 



-^i^Nu = {w'Q'dzQ^^ - (£e,av), (A9) 



(AlO) 
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Adding (jA9p and (jAlOp . the exact relation for (£e)z, 
(|A2[) is obtained. However, this is not the case for the 
relation of the dissipation rate of TKE (jAip . Instead of 
the heat-flux based on the filtered velocity w'Q', (jAip 
contains the heat-flux of the unfiltered velocity as w'9'. 
Substituting w ~ w—djOjdjW (so not correcting for com- 
pressibility effects), the two heat- fluxes are non-trivially 
coupled by 

HTW = WW - d^a^Q'd^w' + a^j{djQ'){djw'). (All) 
Interestingly, the two extra terms on the right-hand 



side correspond to terms typically encountered in the 
transport equation for the turbulent heat-flux [T^ . The 
flrst is normally associated with the molecular diffusive 
transport, and the second with the dissipative cross- 
correlation term, although the sign is opposite here. Av- 
eraging (|Alip over the height and substituting the result 
into (jASp using (|A7[) . gives that (e)^ is given by 



(e), = -^Ra(Nu - + Pg (c? {d,Q'){d,w')) . 

(A12) 
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